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In 1986, Swendsen and Wang proposed a replica Monte Carlo algorithm for spin glasses 
[Phys. Rev. Lett. 57 (1986) 2607]. Two important ingredients are present, (1) the use of 
a collection of systems (replicas) at different of temperatures, but with the same random 
couplings, (2) defining and flipping clusters. Exchange of information between the systems 
is facilitated by fixing the r spin (r = crV 2 ) and flipping the two neighboring systems 
simultaneously. In this talk, we discuss this algorithm and its relationship to replica exchange 
(also known as parallel tempering) and Houdayer's cluster algorithm for spin glasses. We 
review some of the early results obtained using this algorithm. We also present new results 
for the correlation times of replica Monte Carlo dynamics in two and three dimensions and 
compare them with replica exchange. 

§1. Introduction 

The spin glass problem has been studied extensively over the past 30 years. 1 )~ 3 ) 
Monte Carlo simulation has been one of the main tools. However, spin glass models 
are among the most difficult to simulate, due to their extremely slow dynamics at 
low temperatures. Unlike the ferromagnetic models, 4 ) cluster algorithms for spin 
glasses have only limited success. 

In 1986, a cluster algorithm by the name "replica Monte Carlo" was proposed for 
spin glasses. 5 ^ Although several other algorithms were later developed, 6 ^ 9 - 1 replica 
Monte Carlo remains one of the best algorithm for spin glasses. Replica Monte 
Carlo works very well in two dimensions, reducing the correlation time enormously 
comparing to single spin flips. In one dimension, it is similar to the Swendsen- Wang 
cluster algorithm with a constant correlation time independent of system size. In 
three and higher dimensions, it still provides a great improvement over single spin 
flips, although improvement is not as dramatic. As we will show here, replica Monte 
Carlo becomes essentially equivalent to "replica exchange" or "parallel tempering" 
of Hukushima and Nemoto 8 ) in three or higher dimensions. 

§2. Replica Monte Carlo 

Replica Monte Carlo actually contains two ideas, (1) the simultaneous simulation 
of a collection of systems at different temperatures, (2) cluster moves defined through 
the r spin, where r = a^a 1 is the overlap of the systems at nearby temperatures. 
The first idea was later used in parallel tempering and also in a similar proposal by 
Geryer. 10 ) 
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The reason for simulating an array of systems with the same random couplings 
at different temperature Tj is that configurations can be rapidly equilibrated at high 
temperatures and the information can be quickly transferred to low temperatures. 
Exchanges are carried out pairwise for systems at neighboring temperatures. Con- 
sider two replica configurations a 1 and a 2 . The joint distribution of the two spin 
configurations a 1 and a 2 is a product of the two independent Boltzmann distributions 
at temperature T\ and T2, respectively. The combined system has the "Hamiltonian" 

H pair (a\a 2 ) = -J2 {PJijoWj + 2 Jij° 2 i°1) , (2-1) 

where (3 l = l/(kT\) and (3 2 = 1/(^2), with a probability distribution proportional 
to exp(— ff pa i r ((T 1 , cr 2 )) . Clearly, any Monte Carlo moves that keep the above distri- 
bution invariant are valid. Specifically, we can design Monte Carlo moves that satisfy 
detailed balance with respect to the joint distribution. There are many possible such 
choices. 

In our original replica Monte Carlo implementation, 5 )' 11 )' 12 ) we proposed cluster 
moves. The clusters are defined by connected region of the same r-spin, Tj = crjaf. 
Two nearest neighbor sites i and j are considered connected and in the same cluster 
if Ti = Tj. To keep the definition of the r-cluster invariant with respect to the 
move, we must simultaneously flip both systems. More precisely, we rewrite the pair 
Hamiltonian in the form 

flpairOrV) = " £(/?* + P 2 Tirj)Jij<T\a}. (2-2) 

Since r-spins are fixed, the system is a new spin glass with interaction strength pro- 
portional to (3 1 +(3 2 inside a r-cluster and decreased to (3 l —(5 2 between the r-clusters. 
In particular, if (3 l = f3 2 these r clusters are free to flip without costing any energy. 
In such a case, the r clusters also have a good physical meaning - they correspond to 
the Edwards- Anderson order parameter for spin glass. The interactions between the 
clusters are governed by the above Hamiltonian and can be rewritten as an effective 
Hamiltonian between clusters. Let rj a = 1 be the initial spin value for all cluster 
a indicating non- flipping, and rj a = — 1 indicating a flip of the cluster, the effective 
interaction Hamiltonian between clusters is 

H c i s {rj) = - k a,bVaVb, (2-3) 

a,b 

where the effective coupling is from the sum of the couplings along the boundary 
between cluster a and 6, 

k a ,b = (P 1 -P 2 ) £ Jijo}a). (2-4) 

(i£a,j£b) 

Note that the a above should take the value before the cluster flips. We can now sim- 
ulate the above cluster model with any valid Monte Carlo algorithm, e.g., Metropolis 
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single spin flip, or even Swendsen-Wang cluster flip. After the cluster update to vari- 
ables r] a , the original spins are changed according to a\ <— r} a o\ and a 2 <— i] a o~ 2 for 
i £ a. The replica cluster moves must be supplemented by a normal Metropolis 
simulation of a 1 and a 2 to break the conservation of r. The Metropolis step makes 
the algorithm ergodic. 

In our earlier implementation, we made a sequential sweep of all the clusters 
with a heat-bath rate. It is interesting to note that the replica exchange (parallel 
tempering) is equivalent to a trial move where all the r = —1 clusters are flipped, or a 
move where all r = +1 clusters are flipped, followed by a global sign flip. Since the r 
spins are constant with respect to the replica Monte Carlo moves, any move based on 
values of r satisfies detailed balance and is thus valid. Iba gave further discussion for 
the relation between replica Monte Carlo and other extended ensemble methods. 13 ) 
Another connection is with the cluster algorithm of Houdayer. 7 ) In Houdayer's 
implementation, 14 ^ a site is picked at random, the associated r-cluster with that site 
is formed and flipped with probability 1 for systems with same temperature. This 
is similar to the Wolff variation 15 ) of the Swendsen-Wang algorithm. For systems 
between different temperatures, replica exchange was used. 

Another improvement that was suggested in 1986, but not implemented, n ) is to 
introduce a Swendsen-Wang step within the r clusters. In such case, each bond has 
interaction strength (/3 1 + (3 2 )Jijcrjo-j. The presence of a bond is set with probability 
P = 1 — exp(— 2(/3 1 + f3 2 )\Jij\) if the interaction is satisfied, i.e., Jijajaj > 0. The 
bond is absent otherwise. We do not apply this step between the r clusters. This 
further breaks up the r clusters. One technical advantage of having this step is that 
ergodicity can be maintained without evoking a single-spin-flip sweep. Although it 
is a great help at high temperatures, at interesting temperatures for spin glass phase 
transitions, [3 « 1, the value of P « 0.98 is too large. The dynamics is nearly the 
same as the original version. 

As it is obvious, most of the CPU time is devoted to the computation of the 
effective coupling, k a ^. This can be done efficiently at about 0(N) where N is 
the number of spins of the system. This is because that the number of possible 
interactions is bounded by Nd where d is the dimension of the system. We can 
pre-allocate sufficient memory based on cluster size and the total number of clusters. 
In a new implementation of the algorithm, on a 1GHz Itanium 2, each replica Monte 
Carlo sweep (MCS) per spin takes 1.06 fisec (on a 128 x 128 lattice). This should 
be compared with a straightforward Metropolis algorithm, which runs at 0.29 ^usec 
per MCS per spin. Little size dependence is seen on the speed. 

§3. Some Early Results 

The replica Monte Carlo was used to compute two-dimensional spin-glass sus- 
ceptibility 5 ) and low-temperature heat capacity 11 ) for the ±J spin glass model, 
where the coupling Jy takes +J or —J with equal probability. The spin-glass 
susceptibility was analyzed in the form of power-law divergence \ ~ T^, assum- 
ing T c = 0. The exponent 7 5.1 was found. Recent work 7 ) seems to support 
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X ~ exp(2(2 — i~i)(3J) instead, where r\ is the critical exponent associated with corre- 
lation function, g(r) ~ r -7? . 

The low-temperature heat capacity is more interesting. We found 11 - 1 that c ~ 
(3 2 exp(— 2/3 J) from our simulations and argued that the elementary excitations in 
±J Ising spin glass are kink/antikink pairs. This result was disputed, 16 ) but recent 
more extensive calculations with exact algorithms for partition function 17 -* supported 
our result. 

We also used the replica Monte Carlo algorithm with a renormalization group 
analysis. 18 ^ The results show an early confirmation of the absence of a phase 
transition at finite temperature in two dimensions, and clear evidence of a finite- 
temperature spin- glass phase transitions in three and four dimensions. 




Fig. 1. Time correlation function for the order parameter q on a 128 x 128 ±J Ising spin glass 
system, using replica Monte Carlo with /3J distributed from 0.1 to 1.8 in steps of 0.1 (only 0.9 
to 1.8 are plotted). 

§4. Relaxation Times and Comparison 

In order to compute the spin-glass order parameter and susceptibility, two sets 
of replica are used. In each sweep, we do one single-spin-flip sweep for the two sets, 
and replica Monte Carlo between nearby temperatures for both sets, and replica 
Monte Carlo at same temperature between the sets. Swendsen- Wang clusters are 
generated within a r-cluster. The Edwards- Anderson order parameter appears to 
contain the slowest relaxation mode (the energy correlations decay faster). Thus we 
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Fig. 2. Comparison of integrated correlation time on a 128 x 128 lattice for single-spin-flip (tri- 
angles), parallel tempering (squares), and replica Monte Carlo (circles). The K — f3J value is 
distributed from 0.1 to 1.8 in spacing of 0.1. Typically, 10 6 MCS are used. 

compute (for a single realization of random coupling) 



where q = (1/N) \ ^ Tj|, Tj is the overlap spin at the same temperature. 

In figure Q we present correlation functions at some typical temperatures. The 
function has an initial fast decay followed by a slow relaxation at long time. This 
causes a factor of 10 difference between the integrated correlation time, defined by 



and exponential correlation time, f(t) ~ Aexp(—t/r), t — > oo. The initial fast 
decay reflects the quick swap of configurations between neighboring temperatures. 
However, they can be swapped back again, leading to a longer exponential correlation 
times. In any case, the integrated correlation time is relevant to the error in sample 
average. In figure [21 we present T m t for one single random sample of 128 x 128. 
In contrast to single-spin-flip and parallel tempering, remarkably small correlation 
times are observed for the replica Monte Carlo algorithm. 

In figure EJ we compare several algorithms at 32 x 32. Comparing to the results 
in fig. [21 we found that size effect is small. The curve labelled Houdayer is from 
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Fig. 3. Comparison of integrated correlation times on a 32 x 32 lattice for several algorithms. The 
K = (3.J value is distributed from 0.2 to 1.8 in spacing of 0.2. 




Fig. 4. Integrated correlation time on a 12 x 12 x 12 lattice for single-spin-flip, parallel tempering, 
and replica Monte Carlo. The /3J value is distributed from 0.1 to 1.2 in spacing 0.1. 
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a simulation of 2 sets of systems (instead of a large number of sets). A unit time 
step consists of one Metropolis sweep for all the systems for the two sets, single 
cluster moves between the sets, and replica exchanges in a set between neighboring 
temperatures. 

In figure |1J we present the correlation times for a three-dimensional ± J Ising 
spin-glass system of size 12 3 . The replica Monte Carlo includes moves between the 
two sets at same temperature. As one can see, the performance of replica Monte 
Carlo and parallel tempering is nearly the same. This is because in three or higher 
dimensions, there are essentially only two r clusters, one with + r-spin and other — 
r-spin. Flipping either of them is equivalent to exchange of configurations. 

In these correlation time comparisons, we have not fine-tuned the simulation 
parameters (distribution of temperatures, number of sweeps for each type of move, 
and number of sets). Future research on the optimal simulation parameters might 
reveal further improvements in efficiency. 

§5. Conclusion 

We introduced the replica Monte Carlo algorithm and presented new correlation 
time data. These results show that replica Monte Carlo algorithm is extremely 
efficient in two dimensions, and is even much better than parallel tempering. In 
three dimensions, replica Monte Carlo and parallel tempering are of comparable 
efficiency. 
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